Analiza zbioru danych Bodyfat

Słowem wstępu

Dane pochodzą ze strony internatowej Kaggle.com ze zbioru danych Bodyfat.

Dane wymieniają szacunkowe wartości procentowe tkanki tłuszczowej mierzone pod wodą, wagę i pomiary różnych obwodów ciała 252 mężczyzn.

Rows: 252 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (15): Density, BodyFat, Age, Weight, Height, Neck, Chest, Abdomen, Hip, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Density BodyFat Age Weight Height Neck Chest Abdomen Hip Thigh Knee Ankle Biceps Forearm Wrist
1.07 12.3 23 154.25 67.75 36.2 93.1 85.2 94.5 59.0 37.3 21.9 32.0 27.4 17.1
1.09 6.1 22 173.25 72.25 38.5 93.6 83.0 98.7 58.7 37.3 23.4 30.5 28.9 18.2
1.04 25.3 22 154.00 66.25 34.0 95.8 87.9 99.2 59.6 38.9 24.0 28.8 25.2 16.6
1.08 10.4 26 184.75 72.25 37.4 101.8 86.4 101.2 60.1 37.3 22.8 32.4 29.4 18.2
1.03 28.7 24 184.25 71.25 34.4 97.3 100.0 101.9 63.2 42.2 24.0 32.2 27.7 17.7
1.05 20.9 24 210.25 74.75 39.0 104.5 94.4 107.8 66.0 42.0 25.6 35.7 30.6 18.8

Objaśnienie zmiennych zbioru danych:

  • Density - Gęstość wyznaczona na podstawie ważenia pod wodą

  • BodyFat - Procent tkanki tłuszczowej z równania Siri (1956)

  • Age - Wiek (lata)

  • Weight - Waga (w funtach)

  • Height - Wysokość (cale)

  • Neck - Obwód szyi (cm)

  • Chest - Obwód klatki piersiowej (cm)

  • Abdomen - Obwód brzucha (cm)

  • Hip - Obwód bioder (cm)

  • Thigh - Obwód uda (cm)

  • Knee - Obwód kolana (cm)

  • Ankle - Obwód kostki (cm)

  • Biceps - Obwód bicepsa (cm)

  • Forearm - Obwód przedramienia (cm)

  • Wrist - Obwód nadgarstka (cm)

Równanie Siri:

\hspace{0.8cm} Fat(\%) = \frac{495}{density}-450

Przekształcimy wagę z futów na kilogramy i wysokość w calach na cm według standardu środkowo-europejskigo.

Code
dane["Weight2"] <- round(dane$Weight*0.435,2)

dane["Height2"] <- round(dane$Height*2.54,2)

Dodamy kolumnę ID dla łatwiejszej obróbki danych.

Code
ID <- c(1:length(dane$Density))
dane <- cbind("ID" = ID, dane)

Również dodamy nową kolumnę BMI - to wskaźnik masy ciała, który każdy jest w stanie wyznaczyć i zinterpretować w której kategorii znajduje się jego wynik: niedowaga, waga prawidłowa, nadwaga i otyłość. Uważamy, że nowa zmienna BMI będzie dobrym dopełnieniem naszego zbioru danych.

Code
dane["BMI"] <- round(dane$Weight2/(dane$Height2/100)^2)

Dodamy kolumnę category dla zmiennej BMI, która ma 3 stany:

  • za niska
  • w normie
  • za wysoka
Code
dane["Category"] <- as.factor(ifelse(dane$BMI< 20.1, '1',
                              ifelse(dane$BMI< 25.9, '2','3')))

Dodamy też kategorie dla zmiennej BodyFat.

Code
dane["Category2"] <- as.factor(ifelse(dane$BodyFat< 11, '1',
                              ifelse(dane$BodyFat< 24, '2','3')))

Sprawdzimy czy mamy braki danych w poszczególnych kolumnach:

Code
apply(dane, 2, function(x) sum(is.na(x)))
       ID   Density   BodyFat       Age    Weight    Height      Neck     Chest 
        0         0         0         0         0         0         0         0 
  Abdomen       Hip     Thigh      Knee     Ankle    Biceps   Forearm     Wrist 
        0         0         0         0         0         0         0         0 
  Weight2   Height2       BMI  Category Category2 
        0         0         0         0         0 

Nie mamy braków danych w zbiorze. Kolejnym krokiem będzie wyznaczenie podstawowych statystyk (bez kolumny ID):

Code
summary(dane[,-1])
    Density         BodyFat           Age            Weight     
 Min.   :0.995   Min.   : 0.00   Min.   :22.00   Min.   :118.5  
 1st Qu.:1.041   1st Qu.:12.47   1st Qu.:35.75   1st Qu.:159.0  
 Median :1.055   Median :19.20   Median :43.00   Median :176.5  
 Mean   :1.056   Mean   :19.15   Mean   :44.88   Mean   :178.9  
 3rd Qu.:1.070   3rd Qu.:25.30   3rd Qu.:54.00   3rd Qu.:197.0  
 Max.   :1.109   Max.   :47.50   Max.   :81.00   Max.   :363.1  
     Height           Neck           Chest           Abdomen      
 Min.   :29.50   Min.   :31.10   Min.   : 79.30   Min.   : 69.40  
 1st Qu.:68.25   1st Qu.:36.40   1st Qu.: 94.35   1st Qu.: 84.58  
 Median :70.00   Median :38.00   Median : 99.65   Median : 90.95  
 Mean   :70.15   Mean   :37.99   Mean   :100.82   Mean   : 92.56  
 3rd Qu.:72.25   3rd Qu.:39.42   3rd Qu.:105.38   3rd Qu.: 99.33  
 Max.   :77.75   Max.   :51.20   Max.   :136.20   Max.   :148.10  
      Hip            Thigh            Knee           Ankle          Biceps     
 Min.   : 85.0   Min.   :47.20   Min.   :33.00   Min.   :19.1   Min.   :24.80  
 1st Qu.: 95.5   1st Qu.:56.00   1st Qu.:36.98   1st Qu.:22.0   1st Qu.:30.20  
 Median : 99.3   Median :59.00   Median :38.50   Median :22.8   Median :32.05  
 Mean   : 99.9   Mean   :59.41   Mean   :38.59   Mean   :23.1   Mean   :32.27  
 3rd Qu.:103.5   3rd Qu.:62.35   3rd Qu.:39.92   3rd Qu.:24.0   3rd Qu.:34.33  
 Max.   :147.7   Max.   :87.30   Max.   :49.10   Max.   :33.9   Max.   :45.00  
    Forearm          Wrist          Weight2          Height2      
 Min.   :21.00   Min.   :15.80   Min.   : 51.55   Min.   : 74.93  
 1st Qu.:27.30   1st Qu.:17.60   1st Qu.: 69.16   1st Qu.:173.35  
 Median :28.70   Median :18.30   Median : 76.78   Median :177.80  
 Mean   :28.66   Mean   :18.23   Mean   : 77.83   Mean   :178.18  
 3rd Qu.:30.00   3rd Qu.:18.80   3rd Qu.: 85.69   3rd Qu.:183.52  
 Max.   :34.90   Max.   :21.40   Max.   :157.97   Max.   :197.49  
      BMI         Category Category2
 Min.   : 17.00   1: 29    1: 48    
 1st Qu.: 22.00   2:145    2:128    
 Median : 24.00   3: 78    3: 76    
 Mean   : 24.87                     
 3rd Qu.: 26.00                     
 Max.   :159.00                     

Zauważyłyśmy, że wartość minimalna zmiennej BodyFat jest równa 0. Przyjrzymy się tej obserwacji.

Code
dane %>% 
  filter(BodyFat == 0)
   ID Density BodyFat Age Weight Height Neck Chest Abdomen Hip Thigh Knee Ankle
1 182  1.1089       0  40  118.5     68 33.8  79.3    69.4  85  47.2 33.5  20.2
  Biceps Forearm Wrist Weight2 Height2 BMI Category Category2
1   27.7    24.6  16.5   51.55  172.72  17        1         1

Przypuszczamy, że to może być błąd we wprowadzeniu danych.

Code
495/(dane[182,"Density"])-450
[1] -3.611687

Obliczając poziom tłuszczu z równania Siri uzyskałyśmy wartość ujemną, co oznacza, że wzór Siri nie jest optymalny dla każdego przypadku. Ze względu na możliwą do uzyskania wysokość i wagę zostawiamy obserwację w zbiorze danych ze zmienionym poziomem tłuszczu.

Code
dane[182,2] <- 495/(dane[182,"Density"])-450

Przedstawimy podstawowe statystyki danych w bardziej zgrabnej formie:

Code
st_op <- round(apply(dane[,c(1,2,3,8,16,17,18)],2,summary),2)
st_op <- rbind(st_op,St.dev=apply(dane[,c(1,2,3,8,16,17)],2,sd))
colnames(st_op) <- c("Density", "BodyFat","Age" ,"Abdomen", "Weight", "Height", "BMI")
rownames(st_op) <- c('minimum','kwantyl dolny','mediana','średnia','kwantyl górny','maksimum','odchylenie standardowe')
as.data.frame(round(st_op,2)) %>%
kable(caption="Statystyki zbioru")
Statystyki zbioru
Density BodyFat Age Abdomen Weight Height BMI
minimum 1.00 -3.61 0.00 79.30 15.80 51.55 74.93
kwantyl dolny 63.75 1.04 12.47 94.35 17.60 69.16 173.35
mediana 126.50 1.05 19.20 99.65 18.30 76.78 177.80
średnia 126.50 1.04 19.15 100.82 18.23 77.83 178.18
kwantyl górny 189.25 1.07 25.30 105.38 18.80 85.69 183.52
maksimum 252.00 1.10 47.50 136.20 21.40 157.97 197.49
odchylenie standardowe 72.89 0.29 8.37 8.43 0.93 12.78 72.89

Przedstawimy również graficzną prezentację rozkładów :

Code
ggplot(dane, aes(x=Weight2)) + geom_histogram(color="gray47", fill="navajowhite2") + scale_x_continuous(breaks=seq(50, 160, 10)) +
  labs(title = "Histogram wagi [kg]", x = "Waga [kg]", y = "Liczebność") + theme_bw()

Code
ggplot(dane, aes(x=Height2)) + geom_histogram(color="gray47", fill="navajowhite2") + scale_x_continuous(breaks=seq(74, 204, 10)) +
  labs(title = "Histogram wzrostu [cm]", x = "Wzrost [cm]", y = "Liczebność") + theme_bw()

Na wykresach widzimy, że niektóre obserwacje znacznie wychodzą poza tendencję, w dalszym ciągu zobaczymy jak to wpłynie na wyniki.

Ze statystyk opisowych odczytałyśmy wartość maksymalną BMI, która wyniosła 159. Wyświetlimy obserwację posiadającą niemożliwy poziom wskaźnika masy ciała:

Code
dane[which(dane$BMI == 159), c("ID","Weight2", "Height2", "Age")]
   ID Weight2 Height2 Age
42 42   89.17   74.93  44

Jest ona najprawdopodobniej pomyłką we wprowadzaniu danych, przejrzymy się tej obserwacji w dalszych analizach.

Na koniec przedstawimy wykresy, który ilustrują podział wskaźnika masy ciała i poziomu tłuszczu na kategorie: za niska, w normie, za wysoka.

Code
dane[,-1] %>% 
   ggplot(aes(x = Category))+
   geom_bar(color = "gray70", fill = c("wheat1", "wheat2", "wheat3"))+
  scale_fill_manual("legend", values = c("za niska" = "black", "w normie" = "orange", "za wysoka" = "blue"))+
   theme(axis.text.x = c("niska", "w normie", "za wysoka"))+
   labs(title = "Rozkład BMI", x = "kategoria wagi", y = "liczebność")+
   scale_x_discrete(labels = c("za niska","w normie","za wysoka"))+
   guides(fill=guide_legend(title="kategoria"))+
   scale_fill_discrete(breaks=c("1", "2", "3"), labels=c("za niska", "w normie", "za wysoka"))+
   theme_minimal()

Najbardziej liczną jest grupa z prawidłową masą ciała. Tylko 29 osób z badanych 252 mieści się w kategorii BMI poniżej normy.

Code
dane[,-1] %>% 
   ggplot(aes(x = Category2))+
   geom_bar(color = "gray70", fill = c("wheat1", "wheat2", "wheat3"))+
  scale_fill_manual("legend", values = c("za niska" = "black", "w normie" = "orange", "za wysoka" = "blue"))+
   theme(axis.text.x = c("niska", "w normie", "za wysoka"))+
   labs(title = "Rozkład BodyFat", x = "kategoria poziomu tłuszczu", y = "liczebność")+
   scale_x_discrete(labels = c("za niska","w normie","za wysoka"))+
   guides(fill=guide_legend(title="kategoria"))+
   scale_fill_discrete(breaks=c("1", "2", "3"), labels=c("za niska", "w normie", "za wysoka"))+
   theme_minimal()

Rozkład poziomu tłuszczu w grupie badanych mężczyzn zachowuje się podobnie jak rozkład BMI.

Code
dane <- dane[,-21]

Hipotezy

  1. Czy poszczególne obwody: szyja, klatka piersiowa, brzuch, biodra, uda, kolana, kostka, biceps, przedramię i nadgarstek mają wpływ na BMI?

  2. Czy poszczególne miary: obwody klatki piersiowej, brzucha, bioder, uda, kolana, kostki, bicepsa, przedramienia, nadgarstka, wiek, waga i wzrost mają wpływ na poziom tłuszczu?

Funkcje

Code
gauss_markov <- function(mod) { 
  tabelka <- NULL
  wiersz <- NULL
  bp <- bptest(mod)
  gq <- gqtest(mod)
  dw <- dwtest(mod)
  bg <- bgtest(mod, order.by = ~fitted(mod), order = 3)
  sh <- shapiro.test(mod$residuals)
  ks <- ks.test(mod$residuals, "pnorm")
  lil <- nortest::lillie.test(mod$residuals)
  res <- resettest(mod, type = "regressor")
  rai <- raintest(mod)
  har <- harvtest(mod) 
  w <- c(bp$p.value, gq$p.value, dw$p.value, bg$p.value, sh$p.value, ks$p.value, lil$p.value, res$p.value, rai$p.value, har$p.value)
  h0 <- c("jednorodność wariancji", "jednorodność wariancji", "brak autokorelacji reszt","brak autokorelacji do rzędu 3","rozklad reszt jest normalny", "rozklad reszt  jest normalny", "rozklad reszt  jest normalny", "zależność jest liniowa", "zależność jest liniowa", "zależność jest liniowa")
  h1 <- c("brak jednorodności wariancji", "brak jednorodności wariancji", "występuje autokorelacja reszt", "występuje autokorelacja rzędu 3", "rozklad reszt nie jest normalny", "rozklad reszt nie jest normalny", "rozklad reszt nie jest normalny", "brak zależności liniowej", "brak zależności liniowej", "brak zależności liniowej")
  n <- c("Test Breuscha-Pagana", "Test Goldfelda-Quandta", "Test Durbina-Watsona", "Test Breuscha-Godfreya", "Test Shapiro-Wilka", "Test Kolmogorova-Smirnova", "Test Lillieforsa", "Test RESET Ramseya", "Test Rainbow Uttsa", "Test Harveya-Colliera")
  for (i in 1:length(w)) {
    wiersz <- c(n[i], round(w[i], 3), 
                if (w[i]<0.05) {
                  h1[i]} else {h0[i]} )
    tabelka <- rbind(tabelka, wiersz)
  }

  colnames(tabelka) <- c("Nazwa testu", "P-wartość", "Wniosek")
  rownames(tabelka) <- c(1:length(n))
  tabelka %>% 
    kable(caption="Założenia Gaussa Markowa")
}
Code
odstajace <- function(mod, dane, stopien) { 
h <- hatvalues(mod)[hatvalues(mod) > 2*3/nrow(dane)] %>% as.data.frame()
h <- rownames(h)
rsr <- rstandard(mod)[abs(rstandard(mod)) > 2] %>% as.data.frame()
rsr <- rownames(rsr)
rst <- rstudent(mod)[abs(rstudent(mod)) > 3] %>% as.data.frame()
rst <- rownames(rst)
c <- cooks.distance(mod)[cooks.distance(mod) > 4/nrow(dane)] %>% as.data.frame()
c <- rownames(c)
w <- c(h, rsr, rst, c)
b <- unique(w)
counts <- NULL
for (i in 1:length(b)) {
  counts <- c(counts, sum(w==b[i]))
}
counts <- cbind(b, counts)
a <- c("obserwacja", "ile spełnia kryteriów")
counts <- counts[which(counts[,2]>=stopien),] 
colnames(counts) <- a
counts %>%  kable(caption="Obserwacje odstające")
}
Code
odstajace_col <- function(dat){
  a <- c()
for(i in c(7:16, 19)) {
  a <- c(a, dat %>% identify_outliers(colnames(dat)[i]) %>% select(ID))
}
names(a) <- colnames(dat)[c(7:16, 19)]
print(a)
a <- plyr::ldply(a, data.frame)
repetitions <- names(which(table(a[,2])>=5))
print("Obserwacje występujące co najmniej 5 razy: ")
print(repetitions)
}

Hipoteza 1

BMI to inaczej wskaźnik masy ciała - z angielskiego body mass index - jest to współczynnik powstały przez podzielenie masy ciała podanej w kilogramach przez kwadrat wysokości podanej w metrach. Klasyfikacja wskaźnika BMI została opracowana wyłącznie dla dorosłych.

BMI = \frac{masa[kg]}{wysokość^2[m]}

Code
knitr::include_graphics("bmi.jpg")

Model pierwszy

\hat{BMI} = \beta_{0} + \beta_{1} \cdot Neck + \beta_{2} \cdot Chest + \beta_{3} \cdot Abdomen + \beta_{4} \cdot Hip + \beta_{5} \cdot Thigh + \beta_{6} \cdot Knee + \beta_{7} \cdot Ankle + \beta_{8} \cdot Biceps + \beta_{9} \cdot Forearm + \beta{10} \cdot Wrist

Code
model <- lm(BMI~Neck+Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Wrist, data=dane)
summary(model)

Call:
lm(formula = BMI ~ Neck + Chest + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm + Wrist, data = dane)

Residuals:
    Min      1Q  Median      3Q     Max 
 -6.252  -2.111  -0.352   1.163 120.578 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -11.04020   12.45915  -0.886   0.3764  
Neck         -0.56573    0.41818  -1.353   0.1774  
Chest         0.11726    0.17128   0.685   0.4943  
Abdomen       0.06517    0.14352   0.454   0.6501  
Hip           0.48499    0.22593   2.147   0.0328 *
Thigh         0.26534    0.25046   1.059   0.2905  
Knee          0.03295    0.42583   0.077   0.9384  
Ankle        -0.13217    0.40735  -0.324   0.7459  
Biceps       -0.10045    0.31902  -0.315   0.7531  
Forearm       0.03633    0.37193   0.098   0.9223  
Wrist        -1.13466    0.94514  -1.201   0.2311  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.163 on 241 degrees of freedom
Multiple R-squared:  0.241, Adjusted R-squared:  0.2095 
F-statistic: 7.653 on 10 and 241 DF,  p-value: 1.323e-10
Code
summary(dane$BMI-model$fitted.values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -6.252  -2.111  -0.352   0.000   1.163 120.578 

Mediana jest mniejsza od średniej to rozkład reszt jest prawostronnie asymetryczny.

Wykresy diagnostyczne:

Code
autoplot(model)

Założenia Gaussa Markowa

Code
gauss_markov(model)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.039 brak jednorodności wariancji
Test Goldfelda-Quandta 1 jednorodność wariancji
Test Durbina-Watsona 0.769 brak autokorelacji reszt
Test Breuscha-Godfreya 0.794 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0 rozklad reszt nie jest normalny
Test Kolmogorova-Smirnova 0 rozklad reszt nie jest normalny
Test Lillieforsa 0 rozklad reszt nie jest normalny
Test RESET Ramseya 0.034 brak zależności liniowej
Test Rainbow Uttsa 0 brak zależności liniowej
Test Harveya-Colliera 0.128 zależność jest liniowa

Budując model pełny otrzymałyśmy skorygowany R^2 równy 0.21, przyjrzałyśmy się wykresom i zauważyłyśmy, że obserwacja 42 jest nietypowa. Sprawdzimy, czy tak jest na prawdę.

Code
odstajace(model, dane = dane, 2)
Obserwacje odstające
obserwacja ile spełnia kryteriów
42 4
86 2

Z powyższego zestawienia widać, że obserwacja 42 może wpływać na to, że nasz model nie spełnia założeń. Sprawdźmy czy w wyniku usuwania jej ze zbioru, nasz model się poprawi.

Model drugi

Code
dane2 <- dane[-42,]
model2 <- lm(BMI~Neck+Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Wrist, data=dane2)
summary(model2)

Call:
lm(formula = BMI ~ Neck + Chest + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm + Wrist, data = dane2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0303 -0.6120 -0.0121  0.6024  4.0308 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -12.78573    1.66370  -7.685 3.90e-13 ***
Neck          0.07917    0.05612   1.411 0.159594    
Chest         0.13840    0.02287   6.051 5.49e-09 ***
Abdomen       0.11926    0.01917   6.222 2.17e-09 ***
Hip           0.08192    0.03037   2.697 0.007485 ** 
Thigh         0.12382    0.03347   3.700 0.000268 ***
Knee         -0.24602    0.05691  -4.323 2.26e-05 ***
Ankle         0.07697    0.05442   1.414 0.158566    
Biceps        0.06284    0.04262   1.474 0.141674    
Forearm      -0.02033    0.04967  -0.409 0.682587    
Wrist        -0.00730    0.12658  -0.058 0.954061    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.09 on 240 degrees of freedom
Multiple R-squared:  0.9078,    Adjusted R-squared:  0.904 
F-statistic: 236.4 on 10 and 240 DF,  p-value: < 2.2e-16
Code
summary(dane2$BMI-model2$fitted.values)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-3.03026 -0.61203 -0.01207  0.00000  0.60240  4.03083 

Mediana jest bardzo bliska 0, co wskazuję, na brak asymetrii.

Sprawdzimy jak wyglądają wykresy diagnostyczne:

Code
autoplot(model2)

Z wykresów diagnostycznych raczej spodziewamy się niejednorodności wariancji i braku normalności rozkładu reszt. Sprawdzamy to za pomocą testów.

Założenia Gaussa Markowa

Code
gauss_markov(model2)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.001 brak jednorodności wariancji
Test Goldfelda-Quandta 0.753 jednorodność wariancji
Test Durbina-Watsona 0.135 brak autokorelacji reszt
Test Breuscha-Godfreya 0.575 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0.024 rozklad reszt nie jest normalny
Test Kolmogorova-Smirnova 0.814 rozklad reszt jest normalny
Test Lillieforsa 0.15 rozklad reszt jest normalny
Test RESET Ramseya 0 brak zależności liniowej
Test Rainbow Uttsa 0.094 zależność jest liniowa
Test Harveya-Colliera 0.811 zależność jest liniowa

Po sprawdzeniu założeń dochodzimy do wniosku, że jednak większa część z nich na razie nie jest spełniona. Sprawdzimy obserwacje odstające dla poszczególnych zmiennych:

Code
odstajace_col(dane)
$Neck
[1]  39  45 106

$Chest
[1] 39 41

$Abdomen
[1]  39  41 216

$Hip
[1] 35 39 41

$Thigh
[1]  39  41 152 169

$Knee
[1]  39 192 244

$Ankle
[1] 31 39 86

$Biceps
[1] 39

$Forearm
[1]  45 159 175 206 226

$Wrist
[1]  39  41 226 252

$BMI
[1]  39  41  42 216

[1] "Obserwacje występujące co najmniej 5 razy: "
[1] "39" "41"

Zauważyłyśmy, że obserwacje 39, 41 są odstające prawie dla każdej zmiennej. Również obserwacja 39 jest w każdej statystyce rozpatrywanej w funkcji odstajace(). Sprawdzimy jak model będzie wygłądać bez tych obserwacji:

Code
dane3 <- dane[-c(39,41,42),]

Po wyczyszczeniu danych rozdzielamy nasz zbior na testowy i uczący.

Code
ID <- c(1:length(dane$Density))
dane <- cbind("ID" = ID, dane)
set.seed(221)
wektor <- sample(1:249, 180, replace = FALSE)
uczacy <- dane3[wektor,]
testowy <- dane3[-wektor,]

Dalsze analizy będziemy przeprowadzać na zbiorze uczącym.

Model trzeci

Code
model3 <- lm(BMI~Chest+Neck+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm, data=uczacy)
summary(model3)

Call:
lm(formula = BMI ~ Chest + Neck + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm, data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6920 -0.6510  0.0483  0.6423  3.4233 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.163994   1.854714  -4.941 1.85e-06 ***
Chest        0.139161   0.025767   5.401 2.21e-07 ***
Neck        -0.014175   0.059800  -0.237   0.8129    
Abdomen      0.137760   0.022117   6.229 3.57e-09 ***
Hip         -0.002438   0.036103  -0.068   0.9462    
Thigh        0.145343   0.035583   4.085 6.79e-05 ***
Knee        -0.169290   0.071710  -2.361   0.0194 *  
Ankle        0.020266   0.066001   0.307   0.7592    
Biceps       0.042724   0.046566   0.917   0.3602    
Forearm      0.122550   0.061477   1.993   0.0478 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.043 on 170 degrees of freedom
Multiple R-squared:  0.8924,    Adjusted R-squared:  0.8867 
F-statistic: 156.6 on 9 and 170 DF,  p-value: < 2.2e-16

Możemy zauważyć, że prawdziwe BMI odchyla się od linii regresji średnio o około 1 jednostkę.

Code
summary(dane3$BMI-model3$fitted.values)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-11.81280  -3.07204   0.14458   0.06521   2.98281  10.48664 
Code
autoplot(model3)

Wykresy diagnostyczne wyglądają lepiej, niż w poprzednim modelu, spodziewamy się, że większość zalożeń będzie spełnionych.

Założenia Gaussa Markowa

Code
gauss_markov(model3)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.682 jednorodność wariancji
Test Goldfelda-Quandta 0.933 jednorodność wariancji
Test Durbina-Watsona 0.787 brak autokorelacji reszt
Test Breuscha-Godfreya 0.197 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0.867 rozklad reszt jest normalny
Test Kolmogorova-Smirnova 0.98 rozklad reszt jest normalny
Test Lillieforsa 0.843 rozklad reszt jest normalny
Test RESET Ramseya 0.157 zależność jest liniowa
Test Rainbow Uttsa 0.858 zależność jest liniowa
Test Harveya-Colliera 0.908 zależność jest liniowa

W kolejnym kroku po zbudowaniu pełnego modelu przejrzymy się czy da się nasz model ulepszyć, widzimy, że model posiada 9 zmiennych, co jest sporo, chcemy zredukować je ilość za pomocą metody PCA, którą przedstawiamy poniżej:

Najpierw sprawdzamy korelację zmiennych:

Widać, że niektóre zmienne są mocno skorelowane, więc PCA będzie tutaj dobrym wyborem:

Code
model3.pca <- prcomp(uczacy[,-c(1:6,16:20)], scale = T)
summary(model3.pca)
Importance of components:
                          PC1     PC2    PC3    PC4     PC5     PC6     PC7
Standard deviation     2.5164 0.86709 0.7794 0.6490 0.53794 0.49801 0.43170
Proportion of Variance 0.7036 0.08354 0.0675 0.0468 0.03215 0.02756 0.02071
Cumulative Proportion  0.7036 0.78715 0.8547 0.9014 0.93360 0.96116 0.98187
                           PC8     PC9
Standard deviation     0.30138 0.26899
Proportion of Variance 0.01009 0.00804
Cumulative Proportion  0.99196 1.00000
Code
fviz_eig(model3.pca, addlabels = TRUE)

Code
fviz_contrib(model3.pca, choice = "var", axes = 9)

Widzimy, że pierwsze 3 składowych już dałyby nam niemal 86% wyjaśnionej wariancji, co byłoby zgodnie z kryterium.

Tworzymy zmienne PC1,PC2,PC3.

Code
PC1 <- model3.pca$x[,1]
PC2 <- model3.pca$x[,2]
PC3 <- model3.pca$x[,3]
Code
model3_2 <- lm(BMI~PC1+PC2+PC3, data=uczacy)
summary(model3_2)

Call:
lm(formula = BMI ~ PC1 + PC2 + PC3, data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8746 -0.7642  0.0386  0.8813  4.0660 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 24.13333    0.08748 275.868  < 2e-16 ***
PC1         -1.09812    0.03486 -31.500  < 2e-16 ***
PC2         -0.80886    0.10117  -7.995 1.65e-13 ***
PC3         -0.42268    0.11255  -3.756 0.000235 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.174 on 176 degrees of freedom
Multiple R-squared:  0.8588,    Adjusted R-squared:  0.8564 
F-statistic: 356.8 on 3 and 176 DF,  p-value: < 2.2e-16
Code
gauss_markov(model3_2)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.355 jednorodność wariancji
Test Goldfelda-Quandta 0.949 jednorodność wariancji
Test Durbina-Watsona 0.872 brak autokorelacji reszt
Test Breuscha-Godfreya 0.371 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0.334 rozklad reszt jest normalny
Test Kolmogorova-Smirnova 0.362 rozklad reszt jest normalny
Test Lillieforsa 0.433 rozklad reszt jest normalny
Test RESET Ramseya 0.03 brak zależności liniowej
Test Rainbow Uttsa 0.396 zależność jest liniowa
Test Harveya-Colliera 0.587 zależność jest liniowa

R^2 tego modelu to 0.85 Założenia Gaussa-Markowa z różnych testów troszkę się polepszyły.

Model zagnieżdżony

Code
mod_0 <- lm(BMI~1, data = uczacy)
met_tyl <- stats::step(model3, scope = c(mod_0, model3),
                  direction = 'backward', test = 'F', trace = 0)
met_tyl$coefficients
(Intercept)       Chest     Abdomen       Thigh        Knee     Forearm 
 -9.4910775   0.1444103   0.1339823   0.1549096  -0.1639500   0.1378829 
Code
model_zag <- lm(BMI~Chest+Abdomen+Thigh+Knee+Forearm, data=uczacy)
summary(model_zag)

Call:
lm(formula = BMI ~ Chest + Abdomen + Thigh + Knee + Forearm, 
    data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6718 -0.6431  0.0511  0.6350  3.5310 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.49108    1.51928  -6.247 3.11e-09 ***
Chest        0.14441    0.02415   5.981 1.23e-08 ***
Abdomen      0.13398    0.01978   6.774 1.85e-10 ***
Thigh        0.15491    0.02801   5.531 1.15e-07 ***
Knee        -0.16395    0.06147  -2.667  0.00837 ** 
Forearm      0.13788    0.05417   2.545  0.01178 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.033 on 174 degrees of freedom
Multiple R-squared:  0.8918,    Adjusted R-squared:  0.8887 
F-statistic: 286.7 on 5 and 174 DF,  p-value: < 2.2e-16

Z podsumowania można wnioskować, że prawdziwe BMI odchyla się od linii regresji średnio o około 1 jednostkę, tak samo jak w modelu pełnym.

Code
summary(uczacy$BMI-model_zag$fitted.values)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-2.67183 -0.64313  0.05105  0.00000  0.63502  3.53102 

Założenia Gaussa Markowa

Code
gauss_markov(model_zag)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.638 jednorodność wariancji
Test Goldfelda-Quandta 0.918 jednorodność wariancji
Test Durbina-Watsona 0.787 brak autokorelacji reszt
Test Breuscha-Godfreya 0.63 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0.741 rozklad reszt jest normalny
Test Kolmogorova-Smirnova 0.979 rozklad reszt jest normalny
Test Lillieforsa 0.808 rozklad reszt jest normalny
Test RESET Ramseya 0.277 zależność jest liniowa
Test Rainbow Uttsa 0.813 zależność jest liniowa
Test Harveya-Colliera 0.945 zależność jest liniowa
Code
autoplot(model_zag)

Code
anova(model3, model_zag)
Analysis of Variance Table

Model 1: BMI ~ Chest + Neck + Abdomen + Hip + Thigh + Knee + Ankle + Biceps + 
    Forearm
Model 2: BMI ~ Chest + Abdomen + Thigh + Knee + Forearm
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    170 184.79                           
2    174 185.82 -4   -1.0382 0.2388 0.9161

Na podstawie analizy wariancji przyjmujemy, że model prostszy jest lepszy.

Predykcja

Code
pred <- round(predict(model_zag,testowy),1)
fakt <- testowy[,c(1,19)]
Code
cbind(pred, fakt = fakt$BMI) %>% 
ggplot(aes(x = pred, y = fakt))+
  geom_point()+
  geom_smooth()
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Podsumowanie

Postać modelu:

\hat {BMI} = -9.49 + 0.14*Chest + 0.13*Abdomen + 0.15*Thigh - 0.16*Knee + 0.14*Forearm

Nasz model wyjaśnia niemal 89% wariancji. Zmienne wchodzące w skład modelu to obwody: klatka piersiowa, brzuch, uda, kolana i przedramię. Wskazuje to na dużą zależność tych miar ze wskaźnikiem masy ciała. Oznacza to, że dla przeciętnych mężczyzn wskaźnik ten jest wystarczająco dobrą miarą aby określić swój stan fizyczny. Pamiętajmy jednak, że przy dużej objętości tkanki mięśniowej nie jest on miarodajny i raczej wybralibyśmy procentowy wskaźnik poziomu tłuszczu.

Hipoteza 2

Czy poszczególne miary: obwody klatki piersiowej, brzucha, bioder, uda, kolana, kostki, bicepsa, przedramienia, nadgarstka, wiek, waga i wzrost mają wpływ na poziom tłuszczu?

Code
ggplot(dane3,aes(x=BodyFat,y=Weight2)) + geom_point() + theme_bw() + ggtitle("Zależność % tkanki tłuszczowej i wagi") +
  geom_smooth(method='lm') + ggpubr::stat_regline_equation(label.x = 5, label.y = 250)
`geom_smooth()` using formula = 'y ~ x'

Code
ggplot(dane3,aes(x=BodyFat,y=Abdomen)) + geom_point() + theme_bw() + ggtitle("Zależność % tkanki tłuszczowej i obwodu brzucha") +
  geom_smooth(method='lm') + ggpubr::stat_regline_equation(label.x = 5, label.y = 115)
`geom_smooth()` using formula = 'y ~ x'

Model pierwszy

\hat{BodyFat} = \beta_{0} + \beta_{1} \cdot Chest + \beta_{2} \cdot Abdomen + \beta_{3} \cdot Hip + \beta_{4} \cdot Thigh + \beta_{5} \cdot Knee + \beta_{6} \cdot Ankle + \beta_{7} \cdot Biceps + \beta_{8} \cdot Forearm + \beta_{9} \cdot Age + \beta{10} \cdot Weight2 + \beta{11} \cdot Height2+\beta{12} \cdot Wrist

Code
model <- lm(BodyFat~Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Age+Weight2+Height2+Wrist, data=uczacy)
summary(model)

Call:
lm(formula = BodyFat ~ Chest + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm + Age + Weight2 + Height2 + Wrist, 
    data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1699 -3.0621 -0.0413  2.8888  9.0910 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -24.67529   25.91811  -0.952  0.34245    
Chest        -0.05061    0.11979  -0.422  0.67321    
Abdomen       0.81385    0.10564   7.704 1.11e-12 ***
Hip          -0.10127    0.15618  -0.648  0.51761    
Thigh         0.32767    0.16133   2.031  0.04384 *  
Knee          0.32390    0.31165   1.039  0.30017    
Ankle        -0.10546    0.28083  -0.376  0.70774    
Biceps        0.21557    0.18828   1.145  0.25388    
Forearm       0.37294    0.25023   1.490  0.13801    
Age           0.06188    0.03637   1.701  0.09073 .  
Weight2      -0.21411    0.17333  -1.235  0.21847    
Height2      -0.07795    0.08480  -0.919  0.35931    
Wrist        -1.95932    0.61351  -3.194  0.00168 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.097 on 167 degrees of freedom
Multiple R-squared:  0.7645,    Adjusted R-squared:  0.7476 
F-statistic: 45.17 on 12 and 167 DF,  p-value: < 2.2e-16
Code
summary(uczacy$BodyFat-model$fitted.values)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-9.16994 -3.06208 -0.04125  0.00000  2.88881  9.09097 

Założenia Gaussa Markowa

Code
gauss_markov(model)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.866 jednorodność wariancji
Test Goldfelda-Quandta 0.578 jednorodność wariancji
Test Durbina-Watsona 0.987 brak autokorelacji reszt
Test Breuscha-Godfreya 0.859 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0.226 rozklad reszt jest normalny
Test Kolmogorova-Smirnova 0 rozklad reszt nie jest normalny
Test Lillieforsa 0.359 rozklad reszt jest normalny
Test RESET Ramseya 0.352 zależność jest liniowa
Test Rainbow Uttsa 0.845 zależność jest liniowa
Test Harveya-Colliera 0.641 zależność jest liniowa

Model zagniezdżony

Code
mod_0 <- lm(BodyFat~1, data = uczacy)
met_tyl <- stats::step(model, scope = c(mod_0, model),
                  direction = 'backward', test = 'F', trace = 0)
met_tyl$coefficients
 (Intercept)      Abdomen        Thigh      Forearm          Age      Weight2 
-44.07772256   0.83167546   0.41690583   0.51621459   0.07463536  -0.30133107 
       Wrist 
 -1.82469765 
Code
model_zag <- lm(BodyFat~Abdomen+Thigh+Weight2+Wrist, data=uczacy)
summary(model_zag)

Call:
lm(formula = BodyFat ~ Abdomen + Thigh + Weight2 + Wrist, data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.5177 -2.8102 -0.0278  3.1434  9.9275 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -41.85380    9.78803  -4.276 3.12e-05 ***
Abdomen       0.91504    0.06283  14.563  < 2e-16 ***
Thigh         0.34524    0.12056   2.864 0.004701 ** 
Weight2      -0.31348    0.08241  -3.804 0.000196 ***
Wrist        -1.08634    0.52264  -2.079 0.039118 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.143 on 175 degrees of freedom
Multiple R-squared:  0.7476,    Adjusted R-squared:  0.7419 
F-statistic: 129.6 on 4 and 175 DF,  p-value: < 2.2e-16

Możemy zauważyć, że prawdziwy poziom tłuszczu odchyla się od linii regresji średnio o około 4 punkty procentowe.

Code
summary(uczacy$BodyFat-model_zag$fitted.values)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-9.51767 -2.81022 -0.02781  0.00000  3.14340  9.92750 
Code
autoplot(model_zag)

Założenia Gaussa Markowa

Code
gauss_markov(model_zag)
Założenia Gaussa Markowa
Nazwa testu P-wartość Wniosek
Test Breuscha-Pagana 0.811 jednorodność wariancji
Test Goldfelda-Quandta 0.455 jednorodność wariancji
Test Durbina-Watsona 0.974 brak autokorelacji reszt
Test Breuscha-Godfreya 0.168 brak autokorelacji do rzędu 3
Test Shapiro-Wilka 0.334 rozklad reszt jest normalny
Test Kolmogorova-Smirnova 0 rozklad reszt nie jest normalny
Test Lillieforsa 0.25 rozklad reszt jest normalny
Test RESET Ramseya 0.052 zależność jest liniowa
Test Rainbow Uttsa 0.885 zależność jest liniowa
Test Harveya-Colliera 0.718 zależność jest liniowa

Zauważamy, że R dopasowane jest zbliżone w modelu zagnieżdżonym do R dopasowanego w modelu pełnym. Oba modele spełniają wszystkie założenia, jedynie test Kołmogorowa - Smirnova wykazuje brak normalności rozkładu reszt. Wszystkie zmienne wchodzące w skład modelu zagnieżdżonego są istotne statystycznie, więc skłanialibyśmy się ku modelowi zagnieżdżonemu.

Code
anova(model, model_zag)
Analysis of Variance Table

Model 1: BodyFat ~ Chest + Abdomen + Hip + Thigh + Knee + Ankle + Biceps + 
    Forearm + Age + Weight2 + Height2 + Wrist
Model 2: BodyFat ~ Abdomen + Thigh + Weight2 + Wrist
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    167 2803.8                           
2    175 3004.4 -8   -200.65 1.4939 0.1628

Na podstawie analizy wariancji stwierdzamy, że model prostszy jest lepszy.

Predykcja

Code
pred <- round(predict(model_zag,testowy),1)
fakt <- testowy[,c(1,3)]
Code
cbind(pred, fakt = fakt$BodyFat) %>% 
ggplot(aes(x = pred, y = fakt))+
  geom_point()+
  geom_smooth()
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Podsumowanie

Postać modelu:

\hat {BodyFat} = -41.85 + 0.92*Abdomen + 0.35*Thigh - 0.31*Weight2 - 1.09*Wrist

Nasz model wyjśnia ponad 74% wariancji. Miary wchodzące w skład modelu to obwody: brzuch, uda, nadgarstek oraz waga w kilogramach, dzięki czemu każdy (mężczyzna) może łatwo sobie wyznaczyć procentowy poziom tłuszczu i sprawdzić do jakiej kategori należy. Pomaga to w kontrolowaniu kondycji fizycznej, która ma wysoki wpływ na zdrowie człowieka.